Datos satelitales en GRASS GIS

Autor/a

Verónica Andreo

Fecha de publicación

24 de agosto de 2023

Los datos satelitales en general vienen en formato raster >> aplican las mismas reglas.

Los comandos i.* se orientan explícitamente al procesamiento de datos satelitales (aunque algunos puedan usarse para otros datos raster).

Nota

Para más detalles ver el manual Imagery Intro y la wiki Image Processing

Datos para esta sesión

Escenas Landsat 8 (OLI)

Fechas: 14/01/2020 y 02/03/2020
Path/Row: 229/082
CRS: UTM zona 20 N (EPSG:32620)

Descargar las escenas L8 14/01/2020 (979Mb) y L8 02/03/2020 (880Mb) y moverlas a $HOME/gisdata/landsat_data. No descomprimir!

Historia de la mision Landsat

Lanzamientos de satélites Lansat desde 1972

Comparación entre las bandas de todos los satélites Landsat

The Multispectral Scanner System (MSS) aboard Landsats 1–5 had four bands. The Thematic Mapper (TM) aboard Landsats 4 & 5 had seven bands. Landsat 7’s Enhanced Thematic Mapper Plus (ETM+) has 8 bands and Landsats 8 & 9 have 11 bands. The atmospheric transmission values for this graphic were calculated using MODTRAN for a summertime mid-latitude hazy atmosphere (circa 5 km visibility). Fuente: https://landsat.gsfc.nasa.gov/satellites/landsat-9/landsat-9-bands/.

https://www.usgs.gov/landsat-missions

Manos a la obra

Iniciamos GRASS GIS

Iniciamos GRASS GIS en posgar2007_4_cba/PERMANENT

import os

# data directory
homedir = os.path.expanduser('~')

# GRASS GIS database variables
grassdata = os.path.join(homedir, "grassdata")
project = "posgar2007_4_cba"
mapset = "PERMANENT"

# import standard Python packages we need
import sys
import subprocess

# ask GRASS GIS where its Python packages are to be able to run it from the notebook
sys.path.append(
    subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)

Ahora sí, estamos listos para importar los paquetes de GRASS e iniciar una sesión:

# import the GRASS GIS packages we need
import grass.script as gs
import grass.jupyter as gj

# Start the GRASS GIS Session
session = gj.init(grassdata, project, mapset)

Corroboramos la proyección

# check the CRS
gs.read_command("g.proj", flags="p")

Crear un nuevo mapset

Creamos un nuevo mapset llamado landsat8:

# Create a new mapset

# gs.create_environment()
g.mapset -c mapset=landsat8

Listamos los mapsets accesibles

# list all the mapsets in the search path
g.mapsets -p

Listamos los mapas vectoriales disponibles

# list vector maps in all mapsets in the search path
g.list type=vector

Región de interés

Extraemos el radio urbano de Córdoba

# extract Cordoba urban area from `radios_urbanos`
v.extract input=radios_urbanos \
  where="nombre == 'CORDOBA'" \
  output=radio_urbano_cba

Establecemos la región computacional al radio urbano de Córdoba

# set the computational region to the extent of Cordoba urban area
g.region -p vector=radio_urbano_cba

Descargar e importar los datos L8

Instalar la extensión i.landsat:

# install i.landsat toolset
g.extension extension=i.landsat

Buscar escenas de Landsat 8 disponibles

# search for Landsat 8 scenes
i.landsat.download -l settings=$HOME/gisdata/USGS_SETTING.txt \
  dataset=landsat_8_c1 clouds=35 \
  start='2019-10-27' end='2020-03-15'

NO EJECUTAR! Descargar las escenas seleccionadas

# download selected scenes
# i.landsat.download settings=$HOME/gisdata/USGS_SETTING.txt \
#   id=LC82290822020062LGN00,LC82290822020014LGN00 \
#   output=$HOME/gisdata/landsat_data

Imprimir las bandas dentro de la carpeta

# print all landsat bands within landsat_data folder
i.landsat.import -p input=$HOME/gisdata/landsat_data

Imprimir sólo las bandas seleccionadas con un patrón

# print a selection of bands - might be sloooow
i.landsat.import -p \
  input=$HOME/gisdata/landsat_data \
  pattern='B(2|3|4|5|6|8)'

Importar bandas, recortar y reproyectar al vuelo

# import all bands, subset to region and reproject
i.landsat.import -r \
  input=$HOME/gisdata/landsat_data \
  extent=region

Listar bandas importadas y revisar metadatos

# list raster maps
g.list type=raster mapset=.

# check metadata of some imported bands
r.info map=LC08_L1TP_229082_20200114_20200127_01_T1_B4
r.info map=LC08_L1TP_229082_20200114_20200127_01_T1_B8

Opcional: Importar directorio desde la GUI

Select directory
Import settings

Pre-procesamiento de datos satelitales

De número digital (ND) a reflectancia y temperatura

Los datos L8 OLI vienen en 16-bit con rango de datos entre 0 y 65535.
i.landsat.toar convierte ND en reflectancia TOA (y temperatura de brillo) para todos los sensores Landsat. Opcionalmente proporciona reflectancia de superficie (BOA) después de la corrección DOS.
i.atcorr proporciona un método de corrección atmosférica más complejo para gran variedad de sensores (S6).

Definir region computacional a banda de 30m

# set the region to a 30m band
g.region -p raster=LC08_L1TP_229082_20200114_20200127_01_T1_B4

Convertir DN a reflectancia superficial y temperatura - método DOS

# convert from DN to surface reflectance and temperature - requires to uncompress data locally
i.landsat.toar \
  input=LC08_L1TP_229082_20200114_20200127_01_T1_B \
  output=LC08_229082_20200114_toar_B \
  sensor=oli8 \
  metfile=$HOME/gisdata/landsat_data/LC08_L1TP_229082_20200114_20200127_01_T1_MTL.txt \
  method=dos1

Corroborar info antes y después de la conversión para una banda

# list output maps
g.list type=raster mapset=. pattern="*toar*"

# check info before and after for one band
r.info map=LC08_L1TP_229082_20200114_20200127_01_T1_B3
r.info map=LC08_229082_20200114_toar_B3

Banda 10 de L8 con la paleta de colores *kelvin*

Ahora, sigan los mismos pasos para la escena del 02/03/2020. Qué notan de diferente?

Ajuste de color y composiciones RGB

Ajuste de colores para una composición RGB color natural

# enhance the colors
i.colors.enhance \
  red=LC08_229082_20200114_toar_B4 \
  green=LC08_229082_20200114_toar_B3 \
  blue=LC08_229082_20200114_toar_B2 \
  strength=95

Mostrar la combinación RGB - d.rgb

# display RGB
d.mon wx0
d.rgb \
  red=LC08_229082_20200114_toar_B4 \
  green=LC08_229082_20200114_toar_B3 \
  blue=LC08_229082_20200114_toar_B2

Seguir los mismos pasos para una composición falso color 543. Sobre qué bandas debieran realizar el ajuste?

Composiciones color natural 432 y falso color 543

Enmascarado de nubes con banda QA

Landsat 8 proporciona una banda de calidad (QA) con valores enteros de 16 bits que representan las combinaciones de superficie, atmósfera y condiciones del sensor que pueden afectar la utilidad general de un determinado pixel.
La extensión i.landsat.qa reclasifica la banda QA de Landsat 8 de acuerdo a la calidad del pixel.

Nota

Más información sobre la banda QA de L8 en la guía de usuario.

Crear las reglas para identificar las nubes y sombras de nubes

# create a rule set
i.landsat.qa \
  collection=1 \
  cloud_shadow_confidence="Medium,High" \
  cloud_confidence="Medium,High" \
  output=Cloud_Mask_rules.txt

Reclasificar la banda QA en función de las reglas

# reclass the BQA band based on the rule set created
r.reclass \
  input=LC08_L1TP_229082_20200114_20200127_01_T1_BQA \
  output=LC08_229082_20200114_Cloud_Mask \
  rules=Cloud_Mask_rules.txt

Reporte del porcentaje de nubes y sombras

# report % of clouds and shadows
r.report -e map=LC08_229082_20200114_Cloud_Mask units=p

Mostrar el mapa reclasificado

# display reclassified map over RGB
d.mon wx0
d.rgb \
  red=LC08_229082_20200114_toar_B4 \
  green=LC08_229082_20200114_toar_B3 \
  blue=LC08_229082_20200114_toar_B2
d.rast LC08_229082_20200114_Cloud_Mask

Comparar visualmente la cobertura de nubes con la composición RGB 543.

Composición falso color y máscara de nubes

Fusión de datos/Pansharpening

Vamos a usar la banda PAN (15 m) para mejorar la definición de las bandas espectrales de 30 m, por medio de: i.fusion.hpf](https://grass.osgeo.org/grass-stable/manuals/addons/i.fusion.hpf.html), que aplica un método de adición basado en un filtro de paso alto. Otros métodos están implementados en i.pansharpen.

Instalar la extensión i.fusion.hpf

# Install the reqquired addon
g.extension extension=i.fusion.hpf

Cambiar la región a la banda PAN

# Set the region to PAN band (15m)
g.region -p raster=LC08_229082_20200114_toar_B8

Ejecutar la fusión

# Apply the fusion based on high pass filter
i.fusion.hpf -l -c pan=LC08_229082_20200114_toar_B8 \
  msx=`g.list type=raster mapset=. pattern=*_toar_B[1-7] separator=,` \
  suffix=_hpf \
  center=high \
  modulation=max \
  trim=0.0

Listar los mapas resultantes usando un patrón de búsqueda

# list the fused maps
g.list type=raster mapset=. pattern=*_hpf

Visualizar las diferencias con la herramienta mapswipe

# display original and fused maps
g.gui.mapswipe \
  first=LC08_229082_20200114_toar_B5 \
  second=LC08_229082_20200114_toar_B5_hpf

Datos originales 30 m y datos fusionados 15 m

Índices de agua y vegetación

Establecer la máscara de nubes para evitar el cómputo sobre las nubes

# Set the cloud mask to avoid computing over clouds
r.mask raster=LC08_229082_20200114_Cloud_Mask

Calcular el NDVI y establecer la paleta de colores

# Compute NDVI
r.mapcalc \
  expression="LC08_229082_20200114_NDVI = \
  (LC08_229082_20200114_toar_B5_hpf - LC08_229082_20200114_toar_B4_hpf) / \
  (LC08_229082_20200114_toar_B5_hpf + LC08_229082_20200114_toar_B4_hpf) * 1.0"
# Set the color palette
r.colors map=LC08_229082_20200114_NDVI color=ndvi

Calcular NDWI y establecer la paleta de colores

# Compute NDWI
r.mapcalc expression="LC08_229082_20200114_NDWI = \
  (LC08_229082_20200114_toar_B5_hpf - LC08_229082_20200114_toar_B6_hpf) / \
  (LC08_229082_20200114_toar_B5_hpf + LC08_229082_20200114_toar_B6_hpf) * 1.0"
# Set the color palette
r.colors map=LC08_229082_20200114_NDWI color=ndwi

Mostrar los mapas

# display maps in different monitors
d.mon wx0
d.rast map=LC08_229082_20200114_NDVI

d.mon wx1
d.rast map=LC08_229082_20200114_NDWI

NDVI y NDWI a partir de datos Landsat 8

Estimar NDVI y NDWI para la otra escena usando el módulo i.vi

Clasificación No Supervisada

Agrupar las bandas (i.e., hacer un stack): i.group
Generar firmas para n número de clases: i.cluster
Clasificar: i.maxlik

Listar los mapas usando un patrón

# list the bands needed for classification
g.list type=raster mapset=. pattern=*_toar*_hpf

Crear un grupo de imágenes o stack

# add maps to an imagery group for easier management
i.group group=l8 subgroup=l8 \
 input=`g.list type=raster mapset=. pattern=*_toar*_hpf sep=","`

Obtener estadísticos -firmas- para las n clases de interés con una muestra de pixeles

# statistics for unsupervised classification
i.cluster group=l8 subgroup=l8 \
 sig=l8_hpf \
 classes=7 \
 separation=0.6

Realizar la clasificación no supervisada de toda la imagen

# Maximum Likelihood unsupervised classification
i.maxlik group=l8 subgroup=l8 \
 sig=l8_hpf \
 output=l8_hpf_class \
 rej=l8_hpf_rej

Mostrar el mapa clasificado

# display results
d.mon wx0
d.rast map=l8_hpf_class

Información derivada adicional podría obtenerse con los siguientes módulos, entre otros:

medidas de textura: r.texture,
medidas de diversidad: r.diversity,
estadísticas locales con información de contexto: r.neighbors,
transformación tasseled cap: i.tasscap,
etc.

Clasificación en GRASS GIS

Topic classification en los manuales de GRASS GIS
Image classification en la wiki
Ejemplos de clasificación en el curso dictado en Noruega en 2018
Clasificación con Random Forest en la presentación del OpenGeoHub Summer School 2018 en Praga
Detección de cambios con Landsat

Semantic labels

Un concepto bastante nuevo en GRASS GIS son las etiquetas semánticas o semantic labels. Éstas son especialmente relevantes para las imágenes de satélite, ya que nos permiten identificar a qué sensor y banda corresponde una trama determinada. Estas etiquetas son especialmente relevantes a la hora de trabajar con colecciones de imágenes de satélite y también a la hora de clasificar diferentes escenas. Lo veremos más adelante, pero al generar una firma espectral para un determinado conjunto de bandas, puede reutilizarse para clasificar otra escena siempre que las etiquetas semánticas sean las mismas. Cuidado: aunque es posible reutilizar las firmas espectrales para cualquier escena con las mismas bandas, los cambios temporales (estaciones, impacto meteorológico) limitan su aplicabilidad sólo a escenas obtenidas más o menos al mismo tiempo.

--- title: "Datos satelitales en GRASS GIS" author: Verónica Andreo date: today format: html: code-tools: true code-copy: true code-fold: false execute: eval: false cache: false keep-ipynb: true jupyter: python3 --- Los datos satelitales en general vienen en formato raster >> aplican las mismas reglas. Los comandos **i.\*** se orientan explícitamente al procesamiento de datos satelitales (aunque algunos puedan usarse para otros datos raster). :::{.callout-note} Para más detalles ver el manual <a href="https://grass.osgeo.org/grass-stable/manuals/imageryintro.html">Imagery Intro</a> y la wiki <a href="https://grasswiki.osgeo.org/wiki/Image_processing">Image Processing</a> ::: # Datos para esta sesión :::: columns :::{.column width="50%"} Escenas Landsat 8 (OLI) - Fechas: 14/01/2020 y 02/03/2020 - Path/Row: 229/082 - CRS: UTM zona 20 N (EPSG:32620) ::: :::{.column width="50%"} ![](https://landsat.gsfc.nasa.gov/wp-content/uploads/2021/12/ldcm_2012_COL-300x168-1.png){fig-align="right"} ::: :::: Descargar las escenas [L8 14/01/2020 (979Mb)](https://drive.google.com/file/d/1ytQp-xin1FQr_hqtDJRLgK6g4eXwK-WI/view?usp=sharing) y [L8 02/03/2020 (880Mb)](https://drive.google.com/file/d/1Gg8FbhwpIQR-GyYepM4uw_9IOjEnji_N/view?usp=sharing) y moverlas a `$HOME/gisdata/landsat_data`. **No descomprimir!** # Historia de la mision Landsat ![Lanzamientos de satélites Lansat desde 1972](https://landsat.gsfc.nasa.gov/wp-content/uploads/2020-07/Landsat_timeline_20200318_title.gif) ![Comparación entre las bandas de todos los satélites Landsat](https://landsat.gsfc.nasa.gov/wp-content/uploads/2021/12/all_Landsat_bands.png) The Multispectral Scanner System (MSS) aboard Landsats 1–5 had four bands. The Thematic Mapper (TM) aboard Landsats 4 & 5 had seven bands. Landsat 7’s Enhanced Thematic Mapper Plus (ETM+) has 8 bands and Landsats 8 & 9 have 11 bands. The atmospheric transmission values for this graphic were calculated using MODTRAN for a summertime mid-latitude hazy atmosphere (circa 5 km visibility). Fuente: <https://landsat.gsfc.nasa.gov/satellites/landsat-9/landsat-9-bands/>. https://www.usgs.gov/landsat-missions # Manos a la obra ## Iniciamos GRASS GIS Iniciamos GRASS GIS en *posgar2007_4_cba/PERMANENT* ```{python} import os # data directory homedir = os.path.expanduser('~') # GRASS GIS database variables grassdata = os.path.join(homedir, "grassdata") project = "posgar2007_4_cba" mapset = "PERMANENT" ``` ```{python} # import standard Python packages we need import sys import subprocess # ask GRASS GIS where its Python packages are to be able to run it from the notebook sys.path.append( subprocess.check_output(["grass", "--config", "python_path"], text=True).strip() ) ``` Ahora sí, estamos listos para importar los paquetes de GRASS e iniciar una sesión: ```{python} # import the GRASS GIS packages we need import grass.script as gs import grass.jupyter as gj # Start the GRASS GIS Session session = gj.init(grassdata, project, mapset) ``` Corroboramos la proyección ```{python} # check the CRS gs.read_command("g.proj", flags="p") ``` ## Crear un nuevo mapset Creamos un nuevo mapset llamado *landsat8*: ```{python} # Create a new mapset # gs.create_environment() g.mapset -c mapset=landsat8 ``` Listamos los mapsets accesibles ```{python} # list all the mapsets in the search path g.mapsets -p ``` Listamos los mapas vectoriales disponibles ```{python} # list vector maps in all mapsets in the search path g.list type=vector ``` ## Región de interés Extraemos el radio urbano de Córdoba ```{python} # extract Cordoba urban area from `radios_urbanos` v.extract input=radios_urbanos \ where="nombre == 'CORDOBA'" \ output=radio_urbano_cba ``` Establecemos la región computacional al radio urbano de Córdoba ```{python} # set the computational region to the extent of Cordoba urban area g.region -p vector=radio_urbano_cba ``` ## Descargar e importar los datos L8 Instalar la extensión *i.landsat*: ```{python} # install i.landsat toolset g.extension extension=i.landsat ``` Buscar escenas de Landsat 8 disponibles ```{python} # search for Landsat 8 scenes i.landsat.download -l settings=$HOME/gisdata/USGS_SETTING.txt \ dataset=landsat_8_c1 clouds=35 \ start='2019-10-27' end='2020-03-15' ``` NO EJECUTAR! Descargar las escenas seleccionadas ```{python} # download selected scenes # i.landsat.download settings=$HOME/gisdata/USGS_SETTING.txt \ # id=LC82290822020062LGN00,LC82290822020014LGN00 \ # output=$HOME/gisdata/landsat_data ``` Imprimir las bandas dentro de la carpeta ```{python} # print all landsat bands within landsat_data folder i.landsat.import -p input=$HOME/gisdata/landsat_data ``` Imprimir sólo las bandas seleccionadas con un patrón ```{python} # print a selection of bands - might be sloooow i.landsat.import -p \ input=$HOME/gisdata/landsat_data \ pattern='B(2|3|4|5|6|8)' ``` Importar bandas, recortar y reproyectar al vuelo ```{python} # import all bands, subset to region and reproject i.landsat.import -r \ input=$HOME/gisdata/landsat_data \ extent=region ``` Listar bandas importadas y revisar metadatos ```{python} # list raster maps g.list type=raster mapset=. # check metadata of some imported bands r.info map=LC08_L1TP_229082_20200114_20200127_01_T1_B4 r.info map=LC08_L1TP_229082_20200114_20200127_01_T1_B8 ``` #### Opcional: Importar directorio desde la GUI ::: {.panel-tabset} ## Select directory ![](../assets/img/import_directory_1.png){width=50%} ## Import settings ![](../assets/img/import_directory_2.png){width=50%} ::: ## Pre-procesamiento de datos satelitales ![Workflow de pre-procesamiento de datos satelitales](../assets/img/rs_workflow.jpg) ### De número digital (ND) a reflectancia y temperatura - Los datos L8 OLI vienen en 16-bit con rango de datos entre 0 y 65535. - [i.landsat.toar](https://grass.osgeo.org/grass-stable/manuals/i.landsat.toar.html) convierte ND en reflectancia TOA (y temperatura de brillo) para todos los sensores Landsat. Opcionalmente proporciona reflectancia de superficie (BOA) después de la corrección DOS. - [i.atcorr](https://grass.osgeo.org/grass-stable/manuals/i.atcorr.html) proporciona un método de corrección atmosférica más complejo para gran variedad de sensores (S6). Definir region computacional a banda de 30m ```{python} # set the region to a 30m band g.region -p raster=LC08_L1TP_229082_20200114_20200127_01_T1_B4 ``` Convertir DN a reflectancia superficial y temperatura - método DOS ```{python} # convert from DN to surface reflectance and temperature - requires to uncompress data locally i.landsat.toar \ input=LC08_L1TP_229082_20200114_20200127_01_T1_B \ output=LC08_229082_20200114_toar_B \ sensor=oli8 \ metfile=$HOME/gisdata/landsat_data/LC08_L1TP_229082_20200114_20200127_01_T1_MTL.txt \ method=dos1 ``` Corroborar info antes y después de la conversión para una banda ```{python} # list output maps g.list type=raster mapset=. pattern="*toar*" # check info before and after for one band r.info map=LC08_L1TP_229082_20200114_20200127_01_T1_B3 r.info map=LC08_229082_20200114_toar_B3 ``` ![Banda 10 de L8 con la paleta de colores *kelvin*](../assets/img/L8_band10_kelvin.png) :::{.callout--caution title="Tarea"} Ahora, sigan los mismos pasos para la escena del 02/03/2020. Qué notan de diferente? ::: ### Ajuste de color y composiciones RGB Ajuste de colores para una composición RGB color natural ```{python} # enhance the colors i.colors.enhance \ red=LC08_229082_20200114_toar_B4 \ green=LC08_229082_20200114_toar_B3 \ blue=LC08_229082_20200114_toar_B2 \ strength=95 ``` Mostrar la combinación RGB - *d.rgb* ```{python} # display RGB d.mon wx0 d.rgb \ red=LC08_229082_20200114_toar_B4 \ green=LC08_229082_20200114_toar_B3 \ blue=LC08_229082_20200114_toar_B2 ``` :::{.callout--caution title="Tarea"} Seguir los mismos pasos para una composición falso color 543. Sobre qué bandas debieran realizar el ajuste? ::: :::: columns :::{.column width="50%"} ![](../assets/img/L8_composite432.png) ::: :::{.column width="50%"} ![](../assets/img/L8_composite543.png) ::: Composiciones color natural 432 y falso color 543 :::: ### Enmascarado de nubes con banda QA - Landsat 8 proporciona una banda de calidad (QA) con valores enteros de 16 bits que representan las combinaciones de superficie, atmósfera y condiciones del sensor que pueden afectar la utilidad general de un determinado pixel. - La extensión [i.landsat.qa](https://grass.osgeo.org/grass-stable/manuals/addons/i.landsat.qa.html) reclasifica la banda QA de Landsat 8 de acuerdo a la calidad del pixel. :::{.callout-note} Más información sobre la banda QA de L8 en la [guía de usuario](https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/media/files/LSDS-1619_Landsat8-9-Collection2-Level2-Science-Product-Guide-v5.pdf). ::: Crear las reglas para identificar las nubes y sombras de nubes ```{python} # create a rule set i.landsat.qa \ collection=1 \ cloud_shadow_confidence="Medium,High" \ cloud_confidence="Medium,High" \ output=Cloud_Mask_rules.txt ``` Reclasificar la banda QA en función de las reglas ```{python} # reclass the BQA band based on the rule set created r.reclass \ input=LC08_L1TP_229082_20200114_20200127_01_T1_BQA \ output=LC08_229082_20200114_Cloud_Mask \ rules=Cloud_Mask_rules.txt ``` Reporte del porcentaje de nubes y sombras ```{python} # report % of clouds and shadows r.report -e map=LC08_229082_20200114_Cloud_Mask units=p ``` Mostrar el mapa reclasificado ```{python} # display reclassified map over RGB d.mon wx0 d.rgb \ red=LC08_229082_20200114_toar_B4 \ green=LC08_229082_20200114_toar_B3 \ blue=LC08_229082_20200114_toar_B2 d.rast LC08_229082_20200114_Cloud_Mask ``` :::{.callout--caution title="Tarea"} Comparar visualmente la cobertura de nubes con la composición RGB 543. ::: :::: columns :::{.column width="50%"} ![](../assets/img/L8_composite543.png) ::: :::{.column width="50%"} ![](../assets/img/L8_clouds.png) ::: Composición falso color y máscara de nubes :::: ### Fusión de datos/Pansharpening Vamos a usar la banda PAN (15 m) para mejorar la definición de las bandas espectrales de 30 m, por medio de: i.fusion.hpf](https://grass.osgeo.org/grass-stable/manuals/addons/i.fusion.hpf.html), que aplica un método de adición basado en un filtro de paso alto. Otros métodos están implementados en [i.pansharpen](https://grass.osgeo.org/grass-stable/manuals/i.pansharpen.html). Instalar la extensión *i.fusion.hpf* ```{python} # Install the reqquired addon g.extension extension=i.fusion.hpf ``` Cambiar la región a la banda PAN ```{python} # Set the region to PAN band (15m) g.region -p raster=LC08_229082_20200114_toar_B8 ``` Ejecutar la fusión ```{python} # Apply the fusion based on high pass filter i.fusion.hpf -l -c pan=LC08_229082_20200114_toar_B8 \ msx=`g.list type=raster mapset=. pattern=*_toar_B[1-7] separator=,` \ suffix=_hpf \ center=high \ modulation=max \ trim=0.0 ``` Listar los mapas resultantes usando un patrón de búsqueda ```{python} # list the fused maps g.list type=raster mapset=. pattern=*_hpf ``` Visualizar las diferencias con la herramienta mapswipe ```{python} # display original and fused maps g.gui.mapswipe \ first=LC08_229082_20200114_toar_B5 \ second=LC08_229082_20200114_toar_B5_hpf ``` ![Datos originales 30 m y datos fusionados 15 m](../assets/img/L8_mapswipe_hpf.png) ### Índices de agua y vegetación Establecer la máscara de nubes para evitar el cómputo sobre las nubes ```{python} # Set the cloud mask to avoid computing over clouds r.mask raster=LC08_229082_20200114_Cloud_Mask ``` Calcular el NDVI y establecer la paleta de colores ```{python} # Compute NDVI r.mapcalc \ expression="LC08_229082_20200114_NDVI = \ (LC08_229082_20200114_toar_B5_hpf - LC08_229082_20200114_toar_B4_hpf) / \ (LC08_229082_20200114_toar_B5_hpf + LC08_229082_20200114_toar_B4_hpf) * 1.0" # Set the color palette r.colors map=LC08_229082_20200114_NDVI color=ndvi ``` Calcular NDWI y establecer la paleta de colores ```{python} # Compute NDWI r.mapcalc expression="LC08_229082_20200114_NDWI = \ (LC08_229082_20200114_toar_B5_hpf - LC08_229082_20200114_toar_B6_hpf) / \ (LC08_229082_20200114_toar_B5_hpf + LC08_229082_20200114_toar_B6_hpf) * 1.0" # Set the color palette r.colors map=LC08_229082_20200114_NDWI color=ndwi ``` Mostrar los mapas ```{python} # display maps in different monitors d.mon wx0 d.rast map=LC08_229082_20200114_NDVI d.mon wx1 d.rast map=LC08_229082_20200114_NDWI ``` ![](../assets/img/L8_ndvi.png) ![](../assets/img/L8_ndwi.png) NDVI y NDWI a partir de datos Landsat 8 Estimar NDVI y NDWI para la otra escena usando el módulo [i.vi](https://grass.osgeo.org/grass-stable/manuals/i.vi.html) ### Clasificación No Supervisada - Agrupar las bandas (i.e., hacer un stack): [i.group](https://grass.osgeo.org/grass-stable/manuals/i.group.html) - Generar firmas para *n* número de clases: [i.cluster](https://grass.osgeo.org/grass-stable/manuals/i.cluster.html) - Clasificar: [i.maxlik](https://grass.osgeo.org/grass-stable/manuals/i.maxlik.html) Listar los mapas usando un patrón ```{python} # list the bands needed for classification g.list type=raster mapset=. pattern=*_toar*_hpf ``` Crear un grupo de imágenes o *stack* ```{python} # add maps to an imagery group for easier management i.group group=l8 subgroup=l8 \ input=`g.list type=raster mapset=. pattern=*_toar*_hpf sep=","` ``` Obtener estadísticos -firmas- para las *n* clases de interés con una muestra de pixeles ```{python} # statistics for unsupervised classification i.cluster group=l8 subgroup=l8 \ sig=l8_hpf \ classes=7 \ separation=0.6 ``` Realizar la clasificación no supervisada de toda la imagen ```{python} # Maximum Likelihood unsupervised classification i.maxlik group=l8 subgroup=l8 \ sig=l8_hpf \ output=l8_hpf_class \ rej=l8_hpf_rej ``` Mostrar el mapa clasificado ```{python} # display results d.mon wx0 d.rast map=l8_hpf_class ``` ![](../assets/img/L8_unsup_class.png) Información derivada adicional podría obtenerse con los siguientes módulos, entre otros: - medidas de textura: [r.texture](https://grass.osgeo.org/grass-stable/manuals/r.texture.html), - medidas de diversidad: [r.diversity](https://grass.osgeo.org/grass-stable/manuals/addons/r.diversity.html), - estadísticas locales con información de contexto: [r.neighbors](https://grass.osgeo.org/grass-stable/manuals/r.neighbors.html), - transformación tasseled cap: [i.tasscap](https://grass.osgeo.org/grass-stable/manuals/i.tasscap.html), - etc. ### Clasificación en GRASS GIS - [Topic classification](http://grass.osgeo.org/grass-stable/manuals/topic_classification.html) en los manuales de GRASS GIS - [Image classification](http://grasswiki.osgeo.org/wiki/Image_classification) en la wiki - [Ejemplos de clasificación](http://training.gismentors.eu/grass-gis-irsae-winter-course-2018/units/28.html) en el curso dictado en Noruega en 2018 - [Clasificación con Random Forest](https://neteler.gitlab.io/grass-gis-analysis/03_grass-gis_ecad_randomforest/) en la presentación del OpenGeoHub Summer School 2018 en Praga - [Detección de cambios con Landsat](https://veroandreo.gitlab.io/post/jan2021_ilandsat_tutorial/) ## Semantic labels Un concepto bastante nuevo en GRASS GIS son las etiquetas semánticas o *semantic labels*. Éstas son especialmente relevantes para las imágenes de satélite, ya que nos permiten identificar a qué sensor y banda corresponde una trama determinada. Estas etiquetas son especialmente relevantes a la hora de trabajar con colecciones de imágenes de satélite y también a la hora de clasificar diferentes escenas. Lo veremos más adelante, pero al generar una firma espectral para un determinado conjunto de bandas, puede reutilizarse para clasificar otra escena siempre que las etiquetas semánticas sean las mismas. Cuidado: aunque es posible reutilizar las firmas espectrales para cualquier escena con las mismas bandas, los cambios temporales (estaciones, impacto meteorológico) limitan su aplicabilidad sólo a escenas obtenidas más o menos al mismo tiempo.